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We determine the critical temperature of a 3-d homogeneous system of hard-sphere Bosons by 
path-integral Monte Carlo simulations and finite-size scaling. At low densities, we find that the 
critical temperature is increased by the repulsive interactions, in the form of a power law in density 
with exponent 1/3: ATc /To ~ (na^y^^. At high densities the result for liquid helium, namely a 
lower critical temperature than in the non-interacting case, is recovered. We give a microscopic 
explanation for the observed behavior. 
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The observation of Bose-Einstein condensation in 
atomic vapors , at temperatures of a few hundreds of 
nK by evaporative cooling, has revived interest in the 
theoretical investigation of this phenomenon. In this let- 
ter we study bosonic hard-spheres of diameter a confined 
in a cubic box of volume L^. The hard-sphere diameter 
a corresponds to the s-wave scattering length of a real 
inter-atomic potential in a model which is valid at low 
densities and even yields adequate results for relatively 
dense systems such as liquid helium 

We determine the critical temperature Tc of this sys- 
tem for various number densities n and compare it to the 
critical temperature To of the ideal gas j|] . The literature 
provides contradictory results stemming from analytical 
studies, even with regard to the sign of ATc = Tc — To : a 
negative sign is predicted by Hartree-Fock theory Q , and 
a renormalization group calculation higher critical 
temperatures, but with different low density asymptotic 
behaviors, were the result of Ref. §, ATc/Tq - (na^)!/^ 
Ref. 0, ATc/Tq ~ (na^)^^^, and a recent renormaliza- 
tion group calculation which goes beyond the one- 
loop expansion of Ref. | E ] , and which yields ATc /Tq ~ 
(na3)i/6 |. 

We have used a path-integral Monte Carlo method to 
obtain our results. In 3-dimensional Bose systems both 
superfluidity, and Bose-Einstein condensation, are be- 
lieved to occur at the same critical temperature. We de- 
termine this transition temperature by making use of the 
scaling properties of the superfluid fraction ps / p 0,|lll ■ 
From the hypothesis that - close to the bulk transition 
point - thermodynamic quantities in finite systems de- 
pend on temperature only through the ratio of a char- 



acteristic system size, L in our case, and the correlation 
length, ^(t), the following scaling form can be obtained: 



PS 



/p{t) = L-^^'^QiL/m), 



(1) 



valid for relative temperatures t very near the critical 
point: 



Tc 



(2) 



In Eq. (|^) Q{L/S^{t)) is such a function that — for finite 
L — PS / pis analytic in t close tot — 0. ip and —v are the 
bulk critical exponents of the superfluid fraction and the 
correlation length, respectively; their ratio has been mea- 
sured in helium Jl2t a nd calculated with renormalization 
group techniques [[13[; it is consistent with —if/v = — 1. 
At the transition point the correlation length diverges, so 
that Q{L/S^{t)) becomes independent of L. This feature 
allows us to deduce the transition point for the infinite 
system from simulations of finite size samples: curves 
corresponding to various values of the scaled superfluid 
fraction, Lps/p{t), will intersect at i = 0. 

In our path-integral Monte Carlo simulation we mea- 
sure the "stiffness" of the system against twisting the 
phase of the wave function, when a particle is displaced 
across the periodic boundary conditions, by the means of 
the winding number distribution. From that we deduce 
the superfluid fraction . We use the high-temperature 
approximation for the hard-sphere propagator derived in 
Ref. (Our simulations converged too slowly when 

the "image propagator" was employed at a first 
time). 
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For four different particle numbers, N = 27, 64, 125 
and 216, we plot as a function of temperature the scaled 
quantity 



(3) 



in which the dimensionless N^^"^ has replaced the length 
L. 

From Eq. ([|) follows that the four graphs should in- 
tersect at one point which is to be identified as the tran- 
sition point. However, because of the statistical noise in 
our data, an accurate estimate of this point is difficult 
(see Fig. 1). (The same reason prevented us from using 
for an estimator the scaling of the twist free energy as in 
Ref. |jll|l.) According to the scaling hypothesis A^^/'^ps/p 
is analytic close to the transition point. Therefore we fit 
our values of this quantity (for 4 values of N and 10 dif- 
ferent temperatures close to the transition point) to a 
function of the form 

Q(i/^(T))|^^j,^ -Q(0) + giVi/(3-)(T-rc)/rc (4) 

and determine the parameters (5(0), q, v, and Tc- In this 
way we obtain an quantitative estimate for Tc ■ Further- 
more we can check the results by comparing our values 
for Q(0) and the critical exponent v to results obtained 
elsewhere. Assuming that the interacting Bose gas falls 
into the same universality class as the XY model, we ex- 
pect that 



To air/' 



Q{0) 



universal constant 



(5) 



We find values between 0.29 and 0.33 for the system den- 
sities we considered. This is different from the value of 
0.49 ± 0.01 obtained in ||l^. For the critical exponent 
of the correlation length our calculations yield an esti- 
mate —I' = —0.68 ± 0.28, which is in good agreement 
with analytical results jl^ and simulations jl^ for the 
3-d XY model, as well as experimental data which give 
—1/ = —0.67. As a matter of fact, the values for Tc only 
depend slightly on v. results obtained with ly fixed to 
the experimental value differ at most about 0.2% from 
the ones obtained with letting j/ be a free parameter. 
Fig. 2 shows an exemplary fit. 

Fig. 3 shows the results of our calculations, the rela- 
tive change of the critical temperature as a function of 
the density. We distinguish two regimes: at low densities 
the critical temperature for the interacting gas is higher 
than that for the non-interacting system. With increas- 
ing density the critical temperature attains a broad max- 
imum (near na^ w 0.01) before decreasing and coming to 
its minimum value for the highest densities we have used 
in our simulations. For comparison we include in our 
diagram at higher densities experimental and simulation 
results obtained for liquid helium |2^Jll[]. This system is 
well described by an hard-sphere model with an effective 



hard-sphere diameter of a = 2. 2033^4 Superfluid he- 
lium exists only for a limited range of densities before it 
freezes at a density of na^ « 0.24 

We have fitted the values for densities between na^ = 
5 • 10^^ and 5 • lO"'^ against a power law in density 



ATc 
To 



Tc 



Tn 



(6) 



and obtain cq = 0.34 ± 0.06 and 7 = 0.34 ± 0.03. As- 
suming 7 — 1/3, the enhancement of Tc is linear in the 
scattering length a. This exponent is also the result of 
the analytical study in Ref. |Q, which yields, however, a 
value of Co which is about 14 times larger. When mass 
renormalization is included, it is possible that the RG 
calculation in Ref. may also yield this exponent |22) . 

In order to understand our findings on a microscopic 
level, recall the Feynman picture of superfiuidity [ p3| . 
The partition function Z of the system can be written as 
a sum over all states accessible to distinguishable parti- 
cles: 



Trace {S e 

Boltzmann states 



-!3n 



} 



(7) 



where S is the symetrizer with respect to the particle 
labels. We develop the symetrizer S into a sum over 
exchange cycles. At high temperatures the only signifi- 
cant contribution comes from the identity operator (cy- 
cle length 1) and Z reduces to the partition function 
for a classical system. At lower and lower temperatures, 
however, the importance of the contribution of exchange 
cycles of macroscopic length grows. The statistical corre- 
lations extend over longer and longer distances, until - at 
the transition point - a macroscopic structure, identified 
as superfluid component, is established. 

A necessary condition for these exchange cycles to ap- 
pear is that the particles must be close together (sep- 
arated by distances of the order of the thermal wave- 
length At — 27r/i^/(mfcsT)). In an ideal gas, spatial 
density fluctuations are important, and particle clusters 
are likely to appear. This creates regions of lower den- 
sity at other places in the system, which obstruct the 
formation of macroscopic exchange cycles. On the other 
hand, in the case of an moderately dense interacting gas, 
the particles tend to be more homogeneously distributed 
throughout the whole volume of the system. Analogously 
to percolation problems, it is in this case more likely for 
every atom to find a neighbor at a suitable distance, so 
that the exchange cycles can more easily propagate from 
one particle to the next. Superfluidity is hence "easier" 
to achieve, meaning that it can occur at higher temper- 
atures. 

In order to quantify this understanding of the shift of 
the critical temperature, we have determined the two- 
body distribution function g{r). In Fig. 4 we have plot- 
ted g{r) for different densities at the bulk superfluid tran- 
sition temperature for a system of 125 particles. It can be 
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seen easily that at higher densities the "exchange bump" 
of the ideal gas decreases and, because of the hard-sphere 
boundary condition, the distribution tends to become ho- 
mogeneous. Another way to see that the interactions sup- 
press the spatial density fluctuations is to integrate g{r) 
over sub-volumes of the simulation cell which gives the 
fluctuation of the particle number in these sub- volumes 
. We find that the number of particles in sub-spheres 
of radius comparable to \t fluctuates about 30% more 
in the case of an ideal gas than for an interacting gas of 
density 5 • 10~^a~^. 

At higher densities the atoms cannot exchange with 
each other without dragging other particles along. The 
effective mass becomes greater than unity, and the criti- 
cal temperature is decreased as Tc ~ 1 /m* [p|j23| . 

The results obtained in our case of an isotropic sys- 
tem are different from the ones obtained in an harmonic 
external potential |2^j2^ , where the repulsive interac- 
tions decrease the critical temperature at low densities. 
In a trap, interactions give rise to two competing effects. 
The critical temperature is increased by the suppression 
of density fluctuations (even if those are already smaller 
than in the isotropic case because of the confinement). 
On the other hand, the central density is lowered as the 
condensate wave function is broadened. Both terms are 
linear in the scattering length; hence, whether Tc is sup- 
pressed or increased depends on details of the system 
such as the number of particles, the trapping frequency 
U! and the particle mass pTf . 

To conclude, the effect of repulsive interactions is dif- 
ferent for low and high densities: at low densities, the 
"density homogenization" effect prevails so that Tc is 
increased; at high densities, the exchange motion of the 
atoms - which are behind the Bose-Einstein condensation 
phenomenon - are hindered by the hard cores so that the 
critical temperature is decreased. Since, in the experi- 
ments on alkali gases, the value of na^ in the center of 
the trap does not exceed 10~^, the shifts in Tc remain 
small (about 1%), but may be observable in larger traps. 
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Figure captions 

Figure 1: 

Scaled superfluid fraction N^^^ps/ p in function of tem- 
perature T as determined by our simulations for number 
density na^ = 5 • lO'^ and N = 27,64, 125, and 216. Tq 
is the critical temperature of the non-interacting system. 
The temperature where the four lines are crossing is to be 
identified with the critical temperature. A quantitative 
estimate for this temperature is obtained by fitting the 
values for every N with a straight line and determining 
the point of intersection for the fitted curves. 
Figure 2: 

Result of our fit to the data from Fig. 1. t = (T-Tc)/Tc, 
Tc is the critical temperature for the interacting system, 
—V the critical exponent of the correlation length, both 
as obtained by our fit. We find a critical temperature of 
Tc/Ta = 1.057 ±0.002 
Figure 3: 

Critical temperature Tc of an interacting Bose gas ver- 
sus density: a is the hard-sphere diameter, Tq the critical 
temperature of the non-interacting gas. Two different 
scales are used for the vertical axis: for temperatures 
above To the left scale, for values below Tq the right 
scale apply. For comparison, experimental and simu- 
lation results | pO| , pT[ obtained for helium (at densities 



no? = 0.235 and 0.25) are also included. The dashed 
curve at low densities presents a fit of the data points 
between no? ^ h ■ 10"*^ and 5 • 10"^ io I + co{na^)'> , 
yielding 7 = 0.34, the dotted curve is guide to the eye 
(note the change of scale at Tq = To). At zero tempera- 
ture the hard-sphere system freezes at a density of about 
na^ = 0.25 [§. 

Figure 4: 

Two-body distribution function g(r) versus particle dis- 
tance r in units of the thermal wavelength At for different 
densities na^ at the bulk critical temperature. For the 
sake of visibility the single curves are shifted vertically 
by an amout of 0.5. The hard-sphere boundary condi- 
tion changes the short-distance (r — > 0) behaviour dra- 
matically. However, for the "percolation" -like theory of 
superfluidity the behaviour for r « At is more important, 
and here the functions for the ideal and the low density 
gas both feature the appearance of the "exchange bump" 
for r < At: in the non- interacting gas or gases of lower 
density, many particles cluster together at distances less 
than the thermal wavelength. In denser system the inter- 
actions are compensating this statistical attraction and 
the particles fill the system volume homogeneously. This 
allows for superfluidity to appear at higher temperatures. 
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FIG. 2. 
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